clear all
clc
muc=0.033;
C=40;
Lambdac= 2:0.1:15 ;
T=2:0.1:15;
R=zeros(length(T),length(Lambdac));
numbers=zeros(length(T),length(Lambdac));
for k=1:length(T)
for j=1:length(Lambdac)
n=ceil(((Lambdac(j) - muc)*T(k))/(T(k)*C*muc-1));
numbers(k,j)= n;
% Pv=Gaussiandist(n);
Pv=decreasingdist(n);
Pc=Gaussiandist(n);
lb=zeros(n,1);
ub=ones(n,1);
ub(n+1)=inf;
b1=zeros(n+2,1);
lb(n+1)=Lambdac(j) + n*(1/T(k)) ;
for i=1:n
    b1(i)= Lambdac(j)*Pc(i) + (1/T(k)) ;
end
b1(n+1)= -C*n*muc;
b1(n+2)= -muc;
cond = @(x)parameterfun(x,n,Pv,b1);
obj=@(x)x(n+1);
x0 = ones(1,n);
x0(n+1)= Lambdac(j);
[x,fval,exitflag,output,lambda] = fmincon(obj,x0,[],[],[],[],lb,ub,cond);
R(k,j)= x(n+1);
% [cout,ceqout] = parameterfun(x,n,Pv,b1);
% Waitings(j)=max(1/(-cout(1:n)+0.1));
end
end
surf(Lambdac,T,R)

filenames = 'threeD.mat';
save(filenames)

% Lambdac=2:0.1:15;
% T=2:0.1:15;
% for k=1:length(T)
% for j=1:length(Lambdac)
% num(k,j)=ceil(((Lambdac(j) - muc)*T(k))/(T(k)*C*muc-1));
% numc(k,j)=((Lambdac(j) - muc)*T(k))/(T(k)*C*muc-1);
% end
% end
% surf(Lambdac,T,num)
% hold on
% surf(Lambdac,T,numc)